%% Producing Figure 4

close all
warning off

load('no_policy')

%% Generate random shocks

options = optimset;
T = 5000;               % simulation periods
%T = 10000;               % simulation periods
t_drop = 500;

rng('default')
rng(1)                   % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock
%Exi = zeros(T+1);  % stochastic steady state

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
erk  = zeros(T+1,1);
x    = zeros(T+1,npol);
x_nonlog    = zeros(T+1,npol);
x_allnonlog = zeros(T+1,npol);

%% Simulation

zthre = 0;

kv(1)  = stss_k;
rbv(1) = stss_rb;
nv(1)  = stss_n;
lv(1)  = kv(1)/nv(1);
av(1)  = 0;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K_nonlog  = [log(kv(t-1)) rbv(t-1) log(nv(t-1)) xiv(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x(t,:)  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x(t,:)*solved_eta3_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x(t,:)*solved_eta2_rhsee);
            rbv(t)    = exp(x(t,:)*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next - erk(t)) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;

    end
    % End of T period simulation

%% normal times, no runs in the past 12 periods

normal_t = 0;
change_k_normal = zeros(1,1);
change_l_normal = zeros(1,1);
change_n_normal = zeros(1,1);
change_d_normal = zeros(1,1);
p_change_k_normal = zeros(1,1);
p_change_l_normal = zeros(1,1);
p_change_y_normal = zeros(1,1);

for t = t_drop:T

    if sum(runv(t-11:t)) == 0 && pibv(t)>0
        
        normal_t = normal_t + 1;
        change_k_normal(normal_t) = kv(t) - kv(t-1);
        change_l_normal(normal_t) = lv(t) - lv(t-1);
        change_n_normal(normal_t) = nv(t) - nv(t-1);
        change_d_normal(normal_t) = (kv(t)-nv(t)) - (kv(t-1)-nv(t-1));
        
        p_change_k_normal(normal_t ) = (kv(t)/kv(t-1)-1)*100;
        p_change_l_normal(normal_t ) = (lv(t)/lv(t-1)-1)*100;
        p_change_y_normal(normal_t ) = (yv(t)/yv(t-1)-1)*100;
 
    end
    
end

figure('Units','Inches','OuterPosition',[7 0 8 8])
scatter(change_k_normal,change_n_normal,100,'x','LineWidth',1.5)
hold on
scatter(change_k_normal,change_d_normal,'LineWidth',1.5)
xlabel('Change in Asset','Fontsize',30)
ylabel('Change in Equity & Deposit','Fontsize',30)
legend({'Equity','Deposit'},'Fontsize',30)
set(gca,'FontSize',30);
text(0,0,'$\hat{\beta}= 0.98$','Fontsize',30,'Interpreter','latex')   % baseline
grid on

reg_nk  = regstats(change_n_normal,change_k_normal);
reg_dk  = regstats(change_d_normal,change_k_normal);
beta_nk = reg_nk.beta(2)
beta_dk = reg_dk.beta(2)

reg_lk_normal  = regstats(p_change_l_normal,p_change_k_normal);
beta_lk_normal = reg_lk_normal.beta(2)

figure('Units','Inches','OuterPosition',[7 0 8 8])
scatter(p_change_k_normal,p_change_l_normal)
hold on
plot([-6 6],reg_lk_normal.beta(1)+reg_lk_normal.beta(2)*[-6 6],'Linewidth',3)
xlabel('% Change in Asset','Fontsize',30)
ylabel('% Change in Leverage','Fontsize',30)
set(gca,'FontSize',30);
%xlim([-6 6]);
%xticks(-6:2:6)
%ylim([-6 6]);
%yticks(-6:2:6)
text(1,1,'$\hat{\beta} = 0.72$','Fontsize',30,'Interpreter','latex')   % baseline
%text(1,1,'\beta = 0.20','Fontsize',20)   % including negative profit
%text(1,1,'\beta = 0.30','Fontsize',20)   % chi0=0.1
grid on

